#include "SpectrumProcessing.h"

#include <cmath>
#include <iostream>

#ifndef M_PI
    #define M_PI 3.14159265358979323846
#endif

namespace spectClean {

SpectrumProcessing::SpectrumProcessing()
{
}


SpectrumProcessing::~SpectrumProcessing()
{
}

void SpectrumProcessing::dirtySpectrumDFT(const TimeSeries & t, Spectrum & s,
        const double maxFreq, const unsigned overSample)
{
    typedef Spectrum::Complex Complex;

    const unsigned num_times = t.size();
    const double * times = &t.time()[0];
    const double * tAmp = &t.amp()[0];
    const double timeRange = times[num_times-1] - times[0];
    double tMean, tAmpMean;
    t.mean(tMean, tAmpMean);

    std::vector<double> tt(num_times);
    std::vector<double> ta(num_times);
    for (unsigned i = 0; i < num_times; ++i)
    {
        tt[i] = (times[i] - tMean) * 2.0 * M_PI;
        ta[i] = tAmp[i] - tAmpMean;
    }

    const double minFreq = 0.0;
    const double freqRange = maxFreq - minFreq;
    const double freqInc = 1.0 / (timeRange * (double)overSample);
    const unsigned num_freqs = std::ceil(freqRange / freqInc);
    std::cout << " = dirtySpectrumDFT() freqInc = " << freqInc << std::endl;
    std::cout << " = dirtySpectrumDFT() num_freqs = " << num_freqs << std::endl;

    s.resize(2 * num_freqs + 1);
    double * freqs = &s.freq()[0];
    Complex * sAmp = &s.amp()[0];

    // DFT to form spectrum.
#pragma omp parallel for
    for (unsigned j = 0; j < s.size(); ++j)
    {
        freqs[j] = ((double)j - (double)num_freqs) * freqInc;
        sAmp[j] = Complex(0.0, 0.0);

        for (unsigned i = 0; i < num_times; ++i)
        {
            const double arg = -freqs[j] * tt[i];
            sAmp[j] += ta[i] * Complex(std::cos(arg), std::sin(arg));
        }
        sAmp[j] /= num_times;
    }
}




void SpectrumProcessing::windowSpectrumDFT(const TimeSeries & t, Spectrum & s,
        const double maxFreq, const unsigned overSample)
{
    typedef Spectrum::Complex Complex;

    const unsigned num_times = t.size();
    const double * times = &t.time()[0];
    const double * tAmp = &t.amp()[0];
    const double timeRange = times[num_times-1] - times[0];
    double tMean, tAmpMean;
    t.mean(tMean, tAmpMean);

    std::vector<double> tt(num_times);
    std::vector<double> ta(num_times);
    for (unsigned i = 0; i < num_times; ++i)
    {
        tt[i] = (times[i] - tMean) * 2.0 * M_PI;
        ta[i] = tAmp[i] - tAmpMean;
    }

    const double minFreq = 0.0;
    const double freqRange = maxFreq - minFreq;
    const double freqInc = 1.0 / (timeRange * (double)overSample);
    const unsigned num_freqs = std::ceil(freqRange / freqInc);
    std::cout << " = windowSpectrumDFT() freqInc = " << freqInc << std::endl;
    std::cout << " = windowSpectrumDFT() num_freqs = " << num_freqs << std::endl;

    s.resize(2 * num_freqs + 1);
    double * freqs = &s.freq()[0];
    Complex * sAmp = &s.amp()[0];

    // DFT to form spectrum.
#pragma omp parallel for
    for (unsigned j = 0; j < s.size(); ++j)
    {
        freqs[j] = ((double)j - (double)num_freqs) * freqInc;
        sAmp[j] = Complex(0.0, 0.0);

        for (unsigned i = 0; i < num_times; ++i)
        {
            const double arg = -freqs[j] * tt[i];
            sAmp[j] += Complex(std::cos(arg), std::sin(arg));
        }
        sAmp[j] /= num_times;
    }
}



} // namespace spectClean
